module analyticalIS
	use basicmodule
	use setthmod
	implicit none

	integer, parameter	:: nconst=0, nz=nth_one+k+1
	integer				:: aisstage

	contains
	
	function aISgetmsx(x) result(val)
		real	:: x(3),val(3)
		if(stage<3) then
			if(aisstage==0) then
				val=get_msx0t(x,mod(tfcounter,ntailbound)+1)		! search at boundaries first 	
			else
				val=get_msx0t(x,tailflag(1))
			endif
		else
			val=get_msx0t(x,mod(tfcounter,2)+10)		
		endif
	end function
	
	subroutine setanaIS
		use bconf_int
		use nnlpf_int
		use nnlpg_int
		real		:: x(nx0_one),th(nth_one)
		external	:: getminusISratio, aconstr_func, aconstr_grad
		real		:: zlb(nz),zub(nz),z(nz),zguess(nz),zo(nz),r0,r1
		integer, parameter	:: ntries=1000
		integer		:: i,j,nspecial
		
		if(allocated(thlistimp)) deallocate(thlistimp)
		allocate(thlistimp(3,1))
		x=0.7
		tfcounter=1+ntailspecial**2
		thlistimp(:,1)=get_msx0t(x,tailflag(1))
		nimp=1
		tailflag=0
!return		
		zlb=0.001; zub=1-zlb
		call erset(0,0,0)
		do aisstage=0,1
			i=0
			do
				do
					call rnun(zguess)
					if(isvalid(zguess)) exit
				enddo
				call settailflag
	!			call bconf(getminusISratio,0,zlb,zub,x=z,xguess=zguess)
	!			call nnlpf(aconstr_func,3,0,0,zlb,zub,x=z,xguess=zguess,maxitn=50,iprint=0)
				z=0.5
				call nnlpg(aconstr_func,aconstr_grad,nconst,0,0,zlb,zub,x=z,xguess=zguess,maxitn=100,iprint=0)
				i=i+1
				th=aISgetmsx(z(1:3))
				if( getISratio(z)>merge(150,merge(150,200,aisstage==0),stage==3) .and. isvalid(z) .and. all(isfinite(th))) then
!				if( getISratio(z)>merge(150,merge(80,120,aisstage==0),stage==3) .and. isvalid(z) .and. all(isfinite(th))) then
					call mdisp([real(nimp),real(i),getISratio(z)])
!					call mdisp(z)
!					call mdisp([th])
					!block 
					!	real	:: Y(0:k),mv(2)
					!	Y=getYfromZ(th,z)
					!	call mdisp(Y)
					!	call mdisp([getej(1,th),getej(k,th),getej(n0,th)])
					!end block
				
	!				call mdisp(getYfromz(th,z))
					zo=z
			if(th(3)<-0.5) then
				print *,"adjusted the mean"
!						z(1)=min(z(1)+.05,0.999)
!						z(2)=min(z(2)+.05,0.999)
!						th=aISgetmsx(z(1:3))
!				th=[th(1)+0.03/th(2),th(2),th(3)]
			endif
					r0=getISratio(zo)
					nimp=nimp+1
					thlistimp=reshape([thlistimp,th],[3,nimp])
					r1=getISratio(zo)
					print *,"before",r0,"after",r1
					if(r1>400) then
						print *,"adjustment didn't work"
						call mdisp(zo.cvr.z)
						th=aISgetmsx(zo(1:3))
						call mdisp(getYfromZ(th,zo))
						th=aISgetmsx(z(1:3))
						call mdisp(getYfromZ(th,z))
					endif
						
					i=0
					cycle
				endif
				if(i>ntries) exit
			enddo
			print *,"now entering second stage of IS determination"
		enddo
		print *,"done with IS determination"
		x=0
		nspecial=max(nimp/10,1)
		do i=1,ntailspecial
			th=get_msx0t(x,-i)
			do j=1,nspecial
				nimp=nimp+1
				thlistimp=reshape([thlistimp,th],[3,nimp])
			enddo
		enddo
		call mdisp(thlistimp)
		call printtime
	end subroutine
	
	function getYfromz(th,z) result(val)
		use gamin_int
		real	:: th(nth_one),z(nz), val(0:k)
		real	:: eps(0:k),zkom

		eps=z(nx0_one+1:)
		call correps(eps)
		eps(0)=-4+8*eps(0)
		eps(1:k)=-log(eps(1:k))
		call setYfromeps(th,eps,val) 
		
		contains
		
		subroutine correps(eps)
			real	:: eps(0:k),e2,co
			eps=gausscdfinv(eps)
			e2=sum(eps**2)
			co=k+1+5*sqrt(2.0*(k+1))
			if(e2>co) eps=eps/sqrt(e2/co)
			eps=gausscdfvec(eps)
		end subroutine
	end function
		
	function getISratio(z) result(val)
		real	:: z(nz), val
		real	:: x(nx0_one), th(nth_one), Y(0:k), eps(0:k+1)
		integer	:: i
		
		x=z(1:nx0_one)
		th=aISgetmsx(x)
		Y=getYfromz(th,z)
		val=0
		do i=1,nimp
			val=val+getdensimp(thlistimp(:,i),Y)
		enddo
		val=val/nimp
		if(val==0) then
			val=1E10
			return
		endif
		val=getdens0(th,Y)/val
		!if(val/=val) then
		!	Y=getYfromz(th,z)
		!	print *,getdens0(th,Y)
		!	print *,getdens0(th,Y)
		!endif
		
	end function
	
	function isvalid(z) result(val)
		real	:: z(nz), c
		logical	:: val, ierr
		integer	:: i

		do i=1,nconst
			call aconstr_func(z,i,c,ierr)
			if(c<-.01) then
				val=.false.
				return
			endif
		enddo
		val=.true.

	end function
	
	
end module
	
subroutine getminusISratio(n,z,val)
	use analyticalIS
	implicit none
	
	integer	:: n
	real	:: z(n), val
	val=-getISratio(z)
	if(val/=val) val=huge(1.0)
end subroutine
	
subroutine aconstr_func(z,iact,val,ierr)
	use analyticalIS
	implicit none
	integer	:: iact
	real	:: z(nz), val
	logical	:: ierr
	real	:: Y(k),th(nth_one), x(nx0_one)
	
	if(iact==0) then
		val=-getISratio(z)
	else
		val=0
	endif
	ierr=.not. isfinite(val)
	if(ierr) val=1E50
	val=max(-1E10,val)
end subroutine	
	
subroutine aconstr_grad(z,iact,val)
	use analyticalIS
	implicit none
	integer	:: iact
	real	:: z(nz), val(nz)
	logical	:: ierr
	integer	:: i
	real	:: zc(nz), v, vvec(nz)
	real, parameter	:: eps=epsilon(1.0)**(0.33)
!$omp parallel do private(v,zc) num_threads(nz)
	do i=1,nz
		zc=z
		zc(i)=zc(i)+eps
		call aconstr_func(zc,iact,vvec(i),ierr)
		zc(i)=zc(i)-2*eps
		call aconstr_func(zc,iact,v,ierr)
		vvec(i)=vvec(i)-v
		zc(i)=zc(i)+eps
	enddo
	val=vvec/(2*eps)
	if(any(.not. isfinite(val))) val=0
end subroutine	

	